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Abstract 

CN 

^ To understand the computations of our visual system, it is important to understand also the natural 

environment it evolved to interpret. Unfortunately, existing models of the visual environment are either 
^ unrealistic or too complex for mathematical description. Here we describe a naturalistic image model and 

present a mathematical solution for the statistical relationships between the image features and model 
variables. The world described by this model is composed of in<:k;pcndent, opaque, textured objects 
which occlude each other. This simple structure allows us to calculate the joint probability distribution 
,— I of image values sampled at multiple arbitrarily located points, without approximation. This result 

Ch can be converted into probabilistic relationships between observable image features as well as between 

^ the unobservable properties that caused these features, including object boundaries and relative depth. 

Using these results we explain the causes of a wide range of natural scene properties, including highly 
non-gaussian distributions of image features and causal relations between pairs of edges. We discuss the 
implications of this description of natural scenes for the study of vision. 

o 

1 Introduction 

A major goal of vision is to identify physical objects in the world, and their attributes. The relevant 
I— —I sensory evidence — an image — is ambiguous. A visual system must make guesses to interpret this 
^ sensory information, and good guesses should account for the statistics of the input. Consequently, the 
^ statistical structure of natural images has become a subject of fundamental importance for applications 
ranging from computer graphics to neuroscience: Understanding and exploiting natural regularities should 
lead to better visual performance and improved visual representations, whether in image compression or 
p^j in the brain. 

^ Most previous studies of natural scene statistics have characterized natural scenes as linear superposi- 

tions of image features. Principal Components Analysis [1, 2], Independent Components Analysis [3, 4], 
and wavelet transforms [5] each identify related sets of features that when added together can efficiently 
• • reconstitute a natural image. Other methods have improved upon these purely linear, additive descriptions 

. !^ by including multiplicative modulation (e.g. Gaussian scale mixtures [6], hierarchically correlated variances 
[7]). These nonlinear enhancements are useful for representing textures, where common variables like sur- 
^ face properties and illumination intensity and direction naturally co-modulate the contrast of features like 
local orientation. Yet because visual images are not caused by summation but by occlusion, it is important 
to develop models of natural images that are constructed using more accurate nonlinear combinations of 
features [8, 9]. 

We therefore cliosc a simplified model of natural images, colorfully known as the dead leaves model 
[10], for which occlusion is fundamental: The virtual world described by this model is composed of an 
infinitely deep stack of randomly positioned, flat objects ('leaves') that occlude each other (Figure lA- 
B). These objects have attributes of size, shape, texture, and color, independently drawn from specified 
distributions. While the model is only an approximation to our true physical environment, nonetheless it 
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generates images that share many important attributes with natural images, most obviously the ubiquitous 
boundaries between relatively homogeneous regions [11]. It reproduces several known statistical properties 
of natural images, including the spatial power spectrum [11] and bivariate distributions of pixel intensities 
and wavelets coefficients [12]. Despite the demonstrated utility of the model, until now there has been no 
way to calculate higher-order statistical properties of interest except by empirical sampling, which cannot 
provide the insights that exact results can. 

Here we derive an exact solution to the dead leaves model, by calculating joint probability distributions 
explicitly for arbitrary image features. This solution also provides a principled way to relate the features 
in a dead leaves image to the unobserved object attributes that cause these featTires. Since these rela- 
tionships are precisely what we rely upon to see, this result thereby elevates the dead leaves model from 
an interesting approximation of natural images to a valuable tool for modeling perceptual inference and 
neural computation in the visual system. 

To illustrate how this solution helps us understand natural scenes, we apply it to explain the highly 
non-gaussian probability distributions of two important types of image features: wavelet coefficients — 
i.e. the image overlap with localized, oriented filters — and local object boundaries. These features 
are important because they describe stimuli to which neurons in the early visual system are sensitive, 
and l)ccause high-order correlations between them reflect the physical objects and attributes in the visual 
world. Since the functional significance of neural responses to features can depend on the shape of the 
feature distribution [13, 14, 15], it is important to understand why the distributions have their observed 
structure. 

We first look specifically at the marginal, joint, and conditional distributions of wavelet coefficients. 
In natural images, the marginal distributions have heavy tails [16, 17, 6], which we show is due to the 
spatial scale invariance of objects. Joint and conditional distributions of wavelet coefficients have peculiar 
shapes (diamonds, pillows, bowties) that depend on the orientation and distance between the wavelets 
[18, 19, 12]. We show how these distribution shapes arise naturally from occlusion by spatially extended 
objects. Finally, we compute the likelihood that a given pair of local object boundaries comes from the 
same physical contour. When estimated empirically from natural images, this likelihood predicts human 
judgments about contours [20]. Our solution of the dead leaves model recovers the empirical statistics but 
only if one properly accounts for the relative depths at local boundaries, implicating depth cues in simple 
judgments about contours. 

2 Results 

2.1 Solving the dead leaves model 

The pixels of a dead leaves image are fully determined by the properties of objects that are at least partially 
unoccluded. These properties are drawn independently from specified distributions over position, depth, 
size and shape, and texture. Texture can include both mean intensity and (possibly correlated) variations 
about the mean. When we say that we have solved the dead leaves model, we mean that we can calculate 
the joint probabilities of any model variables of interest, whether pixel intensities or object properties. 
This would be straightforward if the image components were related by linear superposition, but is much 
more difficult due to the strong nonlinearity of occlusion. 

The essential property that makes the dead leaves model tractable is that different objects have in- 
dependent attributes. Others have invoked the independence of object properties to derive the two-point 
correlation functions [11] and bivariate intensity probabilities [12] using a recursive argument that accounts 
for the way nearby objects occlude more distant ones. We were able to generalize this calculation from 
two points to an arbitrary collection of N pixels, for which we can now calculate the multivariate joint 
intensity distribution. This distribution can then be transformed into feature probabilities, and related to 
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Figure 1: Example images generated by the dead leaves model. We see layers of objects with random sizes, shapes, 
colors and positions that occlude other objects below. (A) All objects are black or white circles with a relatively 
narrow range of sizes. (B) All objects are textured ellipses with a broad range of sizes drawn from a distribution 
proportional to size^"^, producing approximate scale invariance [11, 12]. Straightforward generalizations allow other 
ensembles of shape and texture. (C) Illustration of an object membership function a. Pixels within a member set of a 
all sample from the same object. Shown is an example dead leaves image with several objects (grey circles) and a set 
of six pixel locations (numbered points). For this configuration, the object membership function is a = {126|3|45}. 

the unobserved object properties. 

If one samples the intensity of a particular dead leaves image at various locations, each pixel value will 
be determined by the texture of whichever object is at the top of the stack at that location. All pixels 
that fall into the same object share its texture, and are thereby correlated; pixels sampling from different 
objects are independent. Thus, if we can specify how the pixels are divided geometrically into objects, 
then we know the complete correlation structure for that image. 

We can mathematically describe the configuration of objects at a given set of N pixels by defining an 
object membership function, a, designating which pixels are 'members' of which objects. (The symbol Q 
is the Hebrew letter mem, chosen to evoke the word membership.) In mathematical language, a is a set 
partition of the N pixels, so it is technically a set of sets: each set corresponds to a different object, and 
it contains the pixel locations at which that object is unobscured by any other objects. For example, one 
might find in a given image that pixels xi, X2 and xg fall into one object, X4 and X5 fall into a different 
object, and X3 is alone in a third object (Figure IC). Then the corresponding object membership function 
can be expressed as a = {{xi,X2,X6}, {X3}, {x4,X5}}, or abbreviated as a = {126|3|45}. 

The object membership function does not contain information directly about the intensities, but only 
about which pixels are correlated. Given a particular object membership a for some selected pixels, the 
probability distribution P(I|a) of image intensities I factorizes into a product over objects: The different 
object textures are independent, and hence so are their respective pixels. In the above example, the 
probability distribution of intensities at those six pixels would be -P(I|a) = P(/i, I2, l6\^)P{l3\^)P{h, -Psja). 
In general, 

|a| 

P(I|a) = n^(I»Ja) (1) 

ra=l 

where |a| is the number of objects, a„ is the set of pixels falling into the nth object of a, and Ia„ is a vector 
of intensities at those pixels. The factors P(Ia„|a) reflect the joint probabilities of intensities in a single, 
textured object. This formulation requires that we specify a texture model to provide these probabilities. 
For concreteness we use a simple gaussian white noise texture superposed on a uniform intensity (Methods) , 
though any other probabilistic texture model could be used instead. Note that the texture model is wholely 
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unrelated to the geometrical aspects of the dead leaves model. 

If the geometric configmation of objects is not known, then the joint distribution of intensities P(l) is 
an average over all possible configurations. The factorized conditional distributions of Equation 1 are then 
combined in the weighted sum 

P(I) = ^P(I|o)P(a) (2) 

a 

This is a mixture distribution in which each mixture component P(I|a) has a distinct correlation structure 
amongst pixels, induced by the different object membership functions. The weighting coefficients are object 
membership probabilities P(Q), i.e. the probability of observing the corresponding memberships over all 
possible dead leaves images with a given shape ensemble. Figure SI shows examples of simple mixture 
distributions. 

The object membership probability P(Q) represents how frequently the N selected pixels are grouped 
into different objects according to 0. We calculate each probability recursively, generalizing an argument 
of [11]. To do so, we must introduce some additional notation. Wc designate as the object membership 
function that remains after removing the nth object. We also define a boolean vector <T(a,n) with A'^ 
components C7i(a,n) = (xj G Q^) that each indicate whether the pixel x, is contained in the nth object of 
a. For instance, (t({126|3|45}, 3) = (0, 0, 0, 1, 1, 0). 

By construction, there is a sequence of objects in any dead leaves image, ordered by depth. Consider 

only the topmost object. There is some probability, which we will denote by Qa-{a,n)^ that this top object 

includes all of the pixels in the set a„, while excluding all the other pixels in ay„. Such an arrangement 

partially satisfies the membership constraint imposed by a. But for this object configuration to contribute 

to P(a), we still need to ensure that the excluded pixels arc also grouped appropriately by objects 'deeper' 

in the image. The probability that deeper objects satisfy these reduced membership constraints is P(ay„). 

Note that this probability is unaffected by whether the deeper objects would have enclosed the pixels in a„: 

Objects at those positions are already occluded by the top object. There is also a probability Qa-{ofl) that 

the top object contains none of the N selected pixels. Given this event, the probability of finding objects 

deeper in the stack that satisfy the membership constraints is just the original factor P(a). Summing 

I qI 

together all possibilities for the top object, we find P(a) = Qa-{a,o)Pi^) + Z^n=i Qcr(a,n)-P('^\n)- Solving for 
P(a) gives the recursion relation 

1 

m = E Q-(a,n)^(a\n) (3) 

Crucially, the image that remains below the top object is yet another dead leaves image, with all the same 
statistical properties as before, so we can calculate P(ay„) by the same formula, recursively. Eventually 
the recursion terminates when there are no pixels left in a, with P(0) = 1. 

This recursive equation applies universally to any dead leaves model with independent, occluding 
objects, regardless of shape. In contrast, the factors Qcr(a,n) depend on the particular shape ensemble and 
the chosen set of pixels. In the Methods section wc derive the general form of these factors for arbitrary 
shapes with smooth boundaries. The Supporting Information (Text SI) provides mathematical details of 
the calculation for the scale-invariant ensemble of elliptical shapes used from here onward. 

The number of possible object membership functions quickly grows large as we consider more pixels. 
The limiting step is the number of possible object membership functions, known as Bell's number B^, 
which unfortunately grows slightly faster than exponentially. In practice this restricts exact calculation to 
around a dozen pixels. Despite this limitation, interesting insights can be gained both by using few pixels 
or few subsets of possible object memberships, and by analyzing the general behavior in various limits. 
For instance, in low-clutter conditions when the maximal distance between pixels u is much smaller than 
the minimum object size r_ (e.g. Figure lA), object membership probabilities P(a) behave as P(a) ~ 
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(u/r_)l"l~^ (Figure S2). Consequently, edges arc rare (lOcdgel = 2), T-junctions are rarer (|aT-junc| = 3), 
and every other feature is rarer still (|a| > 3). In the sections below we use the solution of the dead leaves 
model and relevant approximations to explain complex statistical properties of natural scenes. 

2.2 Feature distributions 

In this section we calculate the joint feature probabilities in specific cases where features are linear functions 
of the pixel intensities, / = FI for some filter matrix F. Wc set the filters of F to be wavelets, local 
derivative operators (edge detectors) with a given orientation and scale. Choosing Haar wavelets, which 
weight intensities by ±1, emphasizes non-gaussianity of feature distributions and thereby establishes a 
more stringent test for the image model [12]. 

It has been previously reported that empirical histograms of different Haar wavelets and wavelet pairs 
in the dead leaves model qualitatively reproduce the marginal and joint distributions in natural scenes 
[12]. Where empirical sampling can, at best, expose these interesting statistical similarities, our analytical 
results let us understand their origins. 

2.2.1 Marginal distributions of wavelet coefficients 

One well-described feature of natural images is that the distribution of spatial derivatives A has heavy 
tails (Figure 2A) well approximated as a generalized Laplace distribution P{A) oc e"'^''^ for an exponent (3 
near 1 [16, 17, 6, 12]. The heavy tails in these distributions cannot be obtained from a standard correlated 
gaussian model, because any projection of a multidimensional gaussian is again gaussian. Higher-order 
statistical structure is required. 

This distribution can be calculated exactly for the dead leaves model by representing the local derivative 
by a simple feature: the intensity difference between nearby points, f = Ii — l2- The resultant feature 
distribution is a mixture of two components, a narrow central peak and a broader tail (Figure 2B,E). While 
this is a more kurtotic distribution than the gaussian texture, it does not closely match natural derivative 
histograms (Figure 2A). 

A simple consideration can account for the discrepancy. In our solution of the dead leaves model, what 
we have described so far as pixels are actually samples at infinitesimal points. In contrast, pixels in natural 
images represent light accumulated over some finite sensor area set by film grain, camera sensor wells, 
or photoreceptor cross-sections. This means that measured pixel values don't directly reflect an intensity 
sampled from an object but instead reflect integrals over unresolved sub-pixel details. When many samples 
are summed over some region Xi, one might naively expect the total /j = &Xi be gaussianly 

distributed. However, the usual central limit theorem does not apply, because of the correlations between 
the variables that is induced by spatially extended objects. These correlations can be segregated by re- 
cxprcssing the total intensity in an image patch as a sum of the mean intensities in visible objects, weighted 
by their visible areas, /j = YlliLi \^n\Jn where is the visible area of the nth object and J„ is the average 
intensity in that area. The summands |On|<^n are approximately independent because they correspond to 
different objects. (They are not strictly independent because the visible areas [Qni are constrained to add 
up to the total area of the patch.) This way of writing F reveals two reasons why the sum docs not converge 
to a normal distribution: the number of summands is a random variable, and the summands themselves 
have long-tailed distributions. 

Scale invariance demands that the areas of homogeneous regions be distributed as a power law 
with exponent 2 [21]. If the mean intensity within an object, Jn, is distributed more narrowly than this, 
then the distribution of visible areas |a„| will dominate the tail behavior of the products |a„| J^. /j is thus a 
sum over a random number (one per visible object) of power-law distributed terms. A generalized central 
limit theorem holds that the distribution of such random sums is a two-sided exponential distribution when 
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summands are power-law distributed with exponents of 2 or higher ([22], Figure 2D).''- Indeed, Figure 2D 
shows nearly straight tails on the log-probability plot. In natural images, areas of homogeneous regions 
are again distributed as power laws, but depending on the particular image the exponents can be slightly 
below 2 [21]. Another generalized central limit theorem shows that under such conditions the distribution 
of the sum has slightly heavier tails [23], as observed (Figure 2A). While these considerations apply directly 
to the average pixel intensities, they pertain equally to intensity differences. 

We can visualize how the heavy tails emerge by plotting feature distributions conditioned on various 
object configurations. When many different objects are visible, the independent object intensities tend 
to average out giving a narrow distribution; when few objects are visible, their areas arc large, so the 
few object intensities are heavily weighted and their distribution is proportionately broad (Figure 2F). 
Components with this spectrum of distribution widths all combine to give the mixture distribution heavy 
tails (Figure 2B-D). 



natural images , dead leaves images 




Figure 2: Log-probability feature distributions log P{f) of spatial derivatives /. (A) Empirically sampled distribution 
of derivatives (depicted graphically, inset) in natural images. (B, C) Feature probabilities calculated exactly for 
dead leaves images, using just one and two samples per image patch respectively (inset). (D) Empirically sampled 
distributions for dead leaves images using a 16x16 grid of samples per patch (inset). (E, F) Mixture components 
P(/|J3) corresponding to panels B and C respectively. 



2.2.2 Joint distributions of wavelet coefficients 

Within natural images, the feature distribution for two orthogonal, colocalizcd wavelets has diamond- 
shaped contours (Figure 3A). Densely sampled dead leaves images reveal the same diamond contours 
(Figure SB); they are already visible when features are represented sparsely (Figure S3A,B). For both 
natural images and dead leaves images, the distinctive non-gaussian structure is most visible for contours 

at large feature amplitude. In this limit, the mixture components with greatest likelihood dominate the 
distribution, and the most likely component at high amplitudes is the one with the greatest variance in the 

^Technically, this theorem requires a geometric distribution for the number of summands. The observed distribution of the 
number of visible objects is not geometric, but it is similarly broad with a width of the same order as its mean, so a similar 
result should hold. 
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given feature direction. The greatest variance for a single Haar wavelet occurs when a boundary between 
two objects aligns with the boundary between oppositely-signed lobes, because that minimizes cancellation 
and maximizes the overlap each object makes with each lobe. However, this arrangement gives a minimal 
variance for the orthogonal wavelet at the same location. As the object boundary rotates (Figure 3C), 
the overlap with one wavelet reduces by exactly the amount that the overlap increases for the orthogonal 
wavelet. Since a mixture component's width is proportional to the overlap, this perfect trade-off gives the 
maximum likelihood contours the observed diamond shape (Methods, Figure 3D).^ 

For neighboring Haar wavelets with the same orientation, there is once again a remarkable similarity 
between pillow-shaped distributions measured for natural images and dead leaves images (Figure 3E,F). 
These too can be explained using simple arguments about the object geometry that dominates at high 
feature amplitudes. Negatively correlated mixture components occur when an object overlaps neighboring 
lobes on neighboring filters (Figure 3G). Positively correlated components occur when an object covers 
the same-signed lobes of both Haar wavelets without cancellation by the intervening lobe of opposite sign. 
This can only happen if a small object occludes that oppositely signed lobe (Figure 3H). Since a very 
limited variety of sizes and positions can accomodate this configuration, the positively correlated mixture 
components have much lower weights. Figure 31 shows mixture components with negative correlations, 
from which emerge the basic 'pillow' shape of the full bivariate feature distributions (Methods) . 

2.2.3 Conditional distributions of wavelet coefiicients 

Wavelet coefficients in natural scenes may be nearly dccorrelated to second order yet still have a strong 
statistical dependency taking the form of a 'bowtie'-shaped distribution of one filter coefficient conditioned 
upon another (Figure 4A) [18, 19]. The dead leaves model reproduces this behavior (Figure 4B), and allows 
us to interpret it as well. 

The distribution of intensities found within an object is narrower than the intensity distribution aver- 
aged over all objects. Consequently, when a wavelet filter lies across an object boundary, it typically yields 
a larger magnitude than the same filter applied wholely within a single object. Since object boundaries 
tend to extend across space, a second filter with different scale or orientation has an elevated probability 
of encountering the same edge. However, as Figs. 4C-D illustrate, the relative sign and magnitude of the 
two feature amplitudes depends on how the object boundary overlaps the second filter. In this symmetric 
example, positive and negative feature amplitudes are equally probable, so the conditional distribution 
-P(/2|/i) broadens with |/i| without any change in the mean (Figure 4E). This explains why the variability 
in one feature amplitude increases with the amplitude of a nearby feature. 

2.3 Shared causes of edges 

A major advantage of using the dead leaves model is that the causes of image features — objects and 
their attributes — are represented explicitly. Our results relate these causes to each other as well as to the 
observable, pixel-based image features. 

In natural images, edge pairs tend to fall tangent to circles passing through both edge locations [24, 25]. 
Geisler et al. [20, 26] augmented such an analysis with global information about physical contours, by 
laboriously hand-segmenting objects within many images of foliage. The likelihood that two edges share 
a physical cause (Figure 5A) — i.e. belong to the same contour — were highly predictive of human 
judgements of whether the edges had a shared cause. 

The dead leaves model can provide a mathematical 'ground truth' for such calculations. First, we 
represent individual edge features by an object membership function that divides four pixels into two pairs 

■^The slight squashing of the diamond shape seen for natural images is a consequence of gravity: in this natural image 
database there are more vertical contours than horizontal ones. Here the dead leaves model does not show this asymmetry 
because it is isotropic by construction. 
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joint log-distribution: example object envelope for 

natural images dead leaves images configurations large features 




Figure 3: Joint feature probabilities logP(/i,/2) for wavelet pairs /i and f-z- For orthogonal, colocalized Haar 
wavelets (inset of A, shifted for visibility), the contours of the empirically sampled bivariate distribution are dia- 
mond shaped for both natural images (A) and dead leaves images (B). At high feature amplitudes, certain object 
configurations have the greatest likelihood and thus dominate the joint distribution. Panel C illustrates one such 
configuration. Colors indicate different objects with unspecified intensities. Dark and light shading show how the 
two wavelets weight the image pixels. (D) Specifying only the object geometry (but not the object intensities) gives 
conditional feature distributions P(/i, f2\^) that are bivariate gaussians with elliptical contours. For the conditional 
distributions that dominate at high feature amplitude, the contours trace a diamond-shaped envelope (thick curve) 
as a function of the relative angles between the object boundary and wavelet orientations (Text 4.3). Parallel, neigh- 
boring wavelets (inset of E) are anticorrelated, with joint probability contours exhibiting a similar 'pillow' shape in 
natural images (E) and dead leaves images (F). Panels G and H illustrate object configurations that dominate at 
large feature amplitudes, colored as in C. (G) If one object covers the opposite-sign lobes of neighboring wavelets 
while others prevent cancellation by the negative lobe, then the conditional feature distribution will have a negative 
correlation. (H) Similarly, if one object covers the same-signed lobes of both features while another object prevents 
cancellation, then the conditional distribution will have a positive correlation. Configurations like this are much 
less probable than those like G, because the middle object must have precisely the right size and position. (I) An 
ensemble of configurations like Panel G produce negatively correlated components (gray ellipses) that vary depending 
on how precisely the objects cover the feature lobes. The positively correlated components (dashed ellipse) caused 
by configurations like Panel H are many times less probable. Discounting the latter gives the mixture distribution 
an overall negative correlation, leaving components that trace out the pillow-shaped envelope (thick curve) seen in 
feature distribution contours (Methods). 
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conditional feature distributions: 
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Figure 4: 'Bowtie' shapes appear in empirically sampled conditional feature distributions -P(/2|/i) for both natural 
images (A) and dead leaves images (B). Horizontal and vertical axes represent the coefficients of two neighboring, 
orthogonal Haar wavelet filters, /i and /2 (inset of A). The grayscale is normalized so black represents and white is 
the maximum probability for a given fi. (C, D) Two equally probable object configurations, colored as in Figure 3C, 
have identical fi but opposite f2- Both features are proportional to the intensity difference between foreground and 
background objects. (E) Conditional feature distribution with only four samples per feature (inset). Traces of the 
limited sampling appear as the faint diagonal bands passing through the origin (highlighted with dotted lines on right 
half). Each distinct band corresponds to a conditional distribution given a different object membership function, 
P(/2|/i, a). Symmetry ensures that there will be no linear correlation between the two features, even as the width of 
P(/2|/i) increases with With features sampled more densely, more such diagonal bands appear, until the bands 
blend together (B). This produces the distinctive bowtie shape in the conditional feature distributions. 



(Figure S4A). Second, we define the conditions under which a pair of edges have the same physical cause. 
Third, for edge pairs with various geometrical relationships (Figure S4B) we plot the likehhood ratio under 
the hypotheses of a shared cause versus different causes (Methods). 

A seemingly natural condition would identify a shared cause when there exists an object that partici- 
pates in both edges. The resultant likelihood ratio always favors a shared cause, for all relative positions 
and orientations of the edge pair (Figure 5B), at odds with reported statistics (Figure 5A) [20, 26]. The 
reason can be seen in Figure 5C: An object could be shared across two edges simply if it is a common 
background for two distinct objects. Thus this definition, only involving object identity on both sides of 
an edge, is inadequate to reproduce the observed edge statistics. 

A more sensible pattern emerges by modifying the definition of common cause to include relative depth, 
assigning 'border ownership' [27] to the local edge. We now define a common cause to exist when a single 
object participates in both edges, and is closer to the viewer than the other objects seen at these edges. An 
example of this configuration is seen in Figure 5E, which agrees with our intuition about a shared cause for 
two edges. Application of this definition requires that the object membership function be augmented to 
include the objects' relative depths, yielding an ordered object membership function. Their probabilities 
can be calculated by a very similar recursion equation as that used for the unordered variant (Methods). 
With this definition. Figure 5D shows that certain edges are more likely to have a common cause, whereas 
other edges are more likely to be independent. The pattern closely resembles results of Geisler et al. 
[20, 26] (Figure 5A). Since those statistics were predictive of human judgments about contour completion 
across occluders, therefore the dead leaves model also qualitatively predicts human inference about such 
ambiguous stimuli. 
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Figure 5: Joint statistics of local edges and global contours. (A) The likelihood ratio that edge pairs in natural 
images are caused by a common object versus by different objects (replotted from [20] with permission). For test 
edges at many distances, directions, and orientations relative to a reference edge (horizontal bar at origin), line 
segments are colored to indicate the likelihood ratio (Methods). The segments are sorted so those indicating high 
likelihoods appear in front. Concentric white rings correspond to unsampled distances. In the dead leaves model, we 
can define the corresponding likelihood in one of two ways. First, a pair of edges could have a 'shared cause' if at 
least one side of each edge samples from the same object. The resultant likelihood is shown in (B) and an example 
of a shared cause is shown in (C). Second, we may add a depth constraint to better describe the existence of a 
shared contour: this shared object must also be on top of the other objects. Using this second definition, panel (D) 
shows the likelihoods and (E) gives an example configuration. These likelihoods reproduce the observations made in 
natural images (A). 
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3 Discussion 



Our study used an occlusion model to explain several distinctive statistical regularities in natural images. 
The model describes images composed of many independent, opaque objects. We solved this image model 
by deriving exact probability distributions that relate arbitrary image features to each other and to the 
depicted objects. By applying and analyzing this solution we were able to account for several curious 
observations about image features, summarized very briefly as follows. We saw that heavy-tailed feature 
distributions are explained by integrating over sub-pixel details with scale-invariant spatial structure (Fig- 
ure 2). The diamond-shaped joint distribution of orthogonal, colocalizcd wavelets occurs because edges 
aligned well with one wavelet must be aligned poorly with an orthogonal wavelet (Figure 3). The pillow- 
shaped joint distribution of parallel wavelets reflects the rarity with which objects can induce positive 
correlation by squeezing precisely into one wavelet lobe (Figure 3). Bowtie-shaped conditional distribu- 
tions arise because extended object boundaries can overlap wavelets with identical amplitudes but opposite 
signs (Figure 4). Finally, accurately computing the likelihood that two edges share a physical cause de- 
pends critically on ascribing relative depth to the edges (Figure 5). The unifying idea is that seemingly 
complex statistics of edge features can be explained by simple geometric configurations of a few opaque 
objects. 

These results were made possible by connecting image features to object configurations through the 
object membership function a. This representation enables probability distributions to be decomposed 
into a mixture of simpler distributions. The existence of a mixture distribution for the dead leaves model 
was first proved in [28, 29]. Here we found an explicit solution for the mixture components that yields 
concrete numbers used in the applications above. Additionally, this solution generalizes to give probabilistic 
relationships among all model variables (Section 4.5), including object texture, size, shape, position, and 
depth. The ability to relate arbitrary image features and many diverse object attributes in a principled 
manner is a substantial advance over previous efforts. 

Although occlusion is a ubiquitous and fundamental attribute of natural scenes, it is not the sole process 
that could cause these effects. However, our results should generalize to other processes that share crucial 
attributes: only one physical cause dominates the image at each point, and separate causes are drawn from 
a scale-invariant size distribution. As one striking example, the cratered lunar surface appears remarkably 
similar to dead leaves images [30]. Even though the causal process is entirely different from occlusion, the 
essential properties are identical: New impacts locally erase traces of previous impacts, and small craters 
are much more common than large ones. Similar principles may approximate other physical processes as 
well, such as those that determine surface composition or some three-dimensional bump textures. The 
results presented here should pertain to feature statistics caused by any such 'exclusion' process. 

3.1 Beyond the dead leaves model 

Despite the dead leaves model's success at reproducing many complex natural statistics, we expect some 
statistical differences also. Indeed, whereas natural scenes appear reasonably gaussian after normalizing 
intensities by the local standard deviation [17, 6], dead leaves images do not have this property. This 
therefore excludes object boundaries as the cause of this property, despite speculations to the contrary 
[31]. By extending the model in various ways, one may hope to capture this and other natural image 
properties and thereby reveal their underlying cause. 

Most real objects have more elaborate shapes than the ellipses used in these calculations. Notably, 
the most common edge configuration seen in natural scenes is consistent with circular [24], elliptical or 
parabolic [26] arcs. This accounts for why the elliptical object ensemble could reproduce statistics of 
images populated by complex, natural objects. Incorporating more complex objects may correct some 
minor discrepancies between the dead leaves model and natural scenes. 

The realism of the dead leaves model could be further improved by adding correlations between model 
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variables. For instance, light sources could be modeled by modulating texture according to position within 
each object. Rudimentary three-dimensional shape could be included using textures to indicate object tilt 
[32] . Perspective could be modeled by covarying size with depth. Binocular disparity could be included by 
generating image pairs in which every object has a positional shift coupled to its depth. Images with such 
improvements could be easily generated, but in some cases a new solution for the enhanced model would 
be required. 

3.2 Toward neural coding of natural scenes 

Some perceptual tasks can be accurately modeled as inference based on simple models of stimulus prob- 
abilities [33, 34, 35, 36, 37]. Human perception of images appears biased toward statistically probable 
features of the dead leaves model. For example, empirical edge statistics predict psychophysical judgments 
about whether two edges have a common cause [20], and the dead leaves model reproduces these statistics. 
Artificial neural networks trained on dead leaves images make systematic interpretation errors that are 
consistent with illusory percepts in humans [38] . Such evidence hints that these percepts might result from 
perceptual inference using probabilities described by the dead leaves model. 

On a more mechanistic level, some electrophysiological recordings of individual neurons in animal cortex 
appear consistent with a probabilistic weighing of sense data [39, 40, 41]. We might speculate that some 
cortical neurons could be tuned to encode feature probabilities. For instance, VI complex cells are excited 
by edges irrespective of polarity and precise location of those edges [42], and are especially sensitive to 
phase alignment caused frequently by object boundaries in natural images [43]. We might therefore wish 
to describe a rudimentary complex cell as encoding the probability that an edge passes through two points 
in its receptive field, irrespective of which side of the edge is brighter. In our formalism, this corresponds to 
an object membership function Oedge = Assuming that objects have gaussian-distribution intensities 

and the image sensors have some additive gaussian noise, the probability of an edge given the intensity 
difference A across space is P(aedgel^) = [l + ^ exp (— ^A^)] , where k and /3 are positive constants that 
depend on the spatial scale, overall image contrast, and sensor noise. This function resembles the contrast- 
energy model of complex cells [44] with a saturating nonlinearity. Thus we might interpret complex cell 
activity as encoding the probability of a local edge in a world of objects. It will be interesting to explore such 
a model more thoroughly, and to see if other neurons have properties that map nicely onto representations 
of still more complex features within the dead leaves model. Since synaptic connections are modified by 
neural correlations, and the occlusion model explains stimulus correlations, therefore the model may also 
help generate predictions about cortical circuitry that has matured in the natural world. 

In vision science, progress has been made by finding stimuli appropriate for the area of study [45]. The 
best stimulus is one that contains a rich repertory of the right kinds of features, while limiting extraneous 
detail. Since the dead leaves model shares many low- and mid-complexity features with the natural 
environment while simplifying some higher-level features, it seems like an especially good stimulus to use 
in experiments that probe the mechanisms of low- and mid-level vision. It strikes a good balance between 
tractability, accuracy, and richness, by isolating two causes of image features which must be disambiguated 
to interpret truly natural scenes: occlusion and texture. The availability of an exact solution for the relevant 
probabilities is a promising new ingredient for experimental and theoretical studies of visual function. 

4 Methods 

4.1 Dead leaves membership probabilities 

Equation 3 expresses the object membership probabilities P(!3) in terms of some geometric factors Qcr(n,n)- These 

factors represent the probability that points Xj G ^„ arc included in one object while the other points G a\„ 
are not, averaged over all object positions and shapes. For convenience, we name these Qo- 'inclusion probabilities'. 
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Note that those quantities involve the geometry of single objects only; the recursion of Equation 3 converts them 
into the multi-object probabilities Pa that characterize the dead leaves model geometry. In this section we show how 
the inclusion probabilities can be calculated for arbitrary objects. 

We begin by specifying a shape through a 'leaf function La{x,p), which is an indicator function over space x 
and shape parameter(s) p. The function can indicate either the inside or the outside of an object centered on the 
origin, depending on the binary variable a € {0, 1}: L„{x,p) equals a when pixel x is inside the object and 1 — a 
when X is outside it (Figure S5A). With this definition, 

N 

<9o-(a,n) (C, P) = n ('^i ~ ^' ^) 

is the inclusion probability that a leaf with shape p and location c includes all sample points x, G a„ and excludes 
all remaining Xj G a\„ (Figure S5B). 

The inclusion probabilities Qcr in Equation 3 are averages over all possible object shapes and positions. Thus we 
are interested in the average of Equation Q„{c, p) over the distribution of leaf positions P(c) and shapes P{p): 



Q^ = j dpP{p)Q^{p) = JJ dpdcP{p)P{c)Q„{c,p) 



We first perform the average over object positions c to obtain Quip), and subsequently calculate the average over 

object shape p. 

In the dead leaves model, objects are distributed with uniform probability across space. For simplicity we also 
assume wraparound boundary conditions and with no loss of generality require that no object is larger than the 
image to avoid self-intersections. (Wc can allow larger objects by choosing a small window into the dead leaves world 
to represent our image; objects may be larger than the window but smaller than the entire model world.) By scaling 
distance so the image has unit area, we have P{c) = 1 and the c-integral of binary- valued Q<t(c, p) gives the inclusion 
probabilities for a given p as the areas of the regions with constant Qcr{c, p). 

Direct integration is not straightforward even for simple object shapes because these regions generally have 
complicated two-dimensional limits and may not even be simply connected. However, using the divergence theorem 
we can transform this area integral into a simpler contour integral that follows object boundaries piecewise. The 
vector field V = has divergence (in c-space) of V • V = 1, so integrating this divergence over the desired region 
gives the enclosed area. The divergence theorem says that this integral equals the flux of V across the region 
boundary: 

Qcr{p) = J dcP(c)0^(c,p) = j V ■Ydc = --ads (4) 

C dC 

where C is the region in c-space where Qcr{c, p) = 1, dC is its boundary, n is the unit normal vector to the boundary, 
and ds is the arclength. The boundary is composed of piecewise smooth segments of the object outline centered 
on the sample points x^ (Figure S5B). We index the relevant segments by m € M, and represent the curves by 
Sm{t) t'^ < t < t'^ iov t between the cusps at which the contour changes direction abruptly. The integral along 
each segment is then 



-r 

2 JtL 



Sm{t) ■ nm{t) ds (5) 



and the complete contour integral is a sum over segments Q^{p) = J2meM 

To average Q^{p) over the shape ensemble P{p) we need to compute J Qrr{p)P{p)dp. Note that the set of 
piecewise smooth segments composing the contour dC may change depending on p, so the p-integral must itself be 
done piecewise. We define an index i specifying the regions in p-space where a given set of segments compose 
the contour. Within R( the integral over p can then be carried out on each summand separately, yielding 



Carrying out this calculation explicitly, not just formally, requires some careful geometry. In the Supporting 
Information we complete these calculations for an ensemble of elliptical objects with an inverse-cube power-law 
distribution of sizes (Text SI). In principle it is also possible to calculate all these probabilities exactly for various 
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other shape ensembles with simple boundaries such as polygons, or compound objects comprising multiple circles. 
Other size ensembles can also be used. The mathematical techniques required to complete the calculations are 
essentially the same. 

For the figures presented in this paper, all objects were ellipses with uniformly distributed eccentricities between 
1 and 4, uniformly distributed orientations, and an inverse-cube size distribution with upper and lower bounds 
r+ = 100 and r_ = 1. For Figs. 2 4, we used high-clutter conditions by setting the pixel spacing to 5r_. For 
Figure 5, to replicate the relatively low-clutter conditions under which the natural image statistics were measured 
empirically [20], we chose the pixel spacing to be r_/5. 

4.2 Intensity and feature distributions 

For simplicity we assume that every object has a constant gaussian-distributed mean intensity and an additive gaus- 

sian white noise tcxtural modulation with variances So and Si. For this texture ensemble, the conditional distribution 
of pixel intensities is oc exp (— |l^C;^^l) , with zero mean and covariance (Co) = Sq X^J^li cri('2, n)aj (Q, n) + 

Si(5jj. In the results shown in this paper, Sq = 1 and Si = 0.01. 

For features specified as linear combinations of intensities by f = FI, the conditional distribution is f (f|l3) oc 
exp (-if^(FCaF^)-if) and the joint probability is the mixture distribution P(f) = -P(!2)-P(f !»)• 

4.3 Averaging over image patches 

Pixels in natural images are integrals of light intensity over a finite solid angle. In the dead leaves model, we can 
approximate these spatial integrals by summing over multiple points within an image patch Xi, defining 

Using the white-noise texture model (Methods 4.2), the total intensity /j over an image patch has a conditional 
distribution which is gaussian with zero mean and variance 

jk n=l 

Here Ca is the covariance matrix of all pixels in image patch Xi conditioned on the object membership function a, 
and is the number of sampled pixels falling into the nth object. Thus the variance increases with the square of 
the sampled area of each object, and is maximized when only one object covers the sampling area. 

A Haar wavelet takes the difference H — Ii — 12 between sums Ii and I2 over two distinct regions (Figure S6A) . 
The corresponding variance does not necessarily increase with the square of each object's sampled area, because 
some of the samples are weighted with opposite signs and thus cancel. The conditional covariance between two Haar 
wavelets Hi and Hj is 

Ch,H,\0 = E - ' \^n' \) ^0 + A^i^'^l (6) 

n=l 

where is the number of samples in region k of wavelet i which fall into the nth object (Figure S6A), and Nij is 
the number of samples shared by wavelets Hi and Hj. 

In Figure 3B,D, the diamond-shaped contours emerge as a consequence of Equation 6. Instead of the Haar 
wavelets with square support shown in that figure, it is simpler to understand the case with circular support (Figure 
S6B), though the result is the same. The maximum amplitude features occur when a single object boundary passes 
through the center of the wavelet at an angle 6. The covariance of the mixture distribution conditioned on this object 
configuration is 

r _oAr2^ 20(7r-20)\ 

where N is the number of samples in each Haar wavelet. For large N this covariance matrix is nearly singular, with 
almost unity correlation coefficient between the variations along Hi and H2 ■ Contours of the corresponding bivariate 
gaussian have maximum extent at feature amplitudes proportional to (±^, ±{^ — 6)). The envelope of these contours 
produces the diamond shown in Figure 3B,D. 
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In Figure 3, two neighboring, parallel Haar wavelets have a joint distribution with a distinctive 'pillow' shape. 
The dominant contributions at high feature amplitudes involve three objects as depicted in Figure 3G, one covering 
the left edge of the wavelet, a second one covering the right edge, and a third covering the gap between them. We 
can approximate this arrangement with a one-dimensional version, considering only the horizontal extent of objects 
(Figure S6C,D). If we denote how much the leftmost and rightmost objects overlap the wavelets by di and dr, then 
the covariance of the mixture distribution is 

where A, = min(rfi,l — di) and the width of each lobe of the Haar wavelet is 1. These components all have a 
correlation coefficient of nearly ±1/2 but have different variances. By changing di and rf,. for the configuration shown 
in Figure S6C we obtain conditional distributions with the ensemble of contours seen in Figure 31. Their envelope 
produces the 'pillow' shape (Figure 3). 



4.4 Shared causes of edges 

To define oriented edges, we select four pixels arranged in a rectangle, and select only those object membership func- 
tions that bisect these four pixels into two pairs. Note that a range of object boundaries can produce such a separation. 
Giving the rectangle an aspect ratio 2.75 constrains edges to an allowed range of orientations 2tan~^ (1/2.75) = 40° 
(Figure S4A) that matches the orientation bandwidth of used in [20]. Pairs of edges are described by two such 
bisected four-pixel clusters (Figure S4B). This definition of edge pairs restricts these eight pixels to have one of only 
seven possible object membership functions (Table lA). In one of these configurations, every pixel pair is a member 
of a different object: a = {12|34|56|78}. In the remaining configurations, at least two pairs are members of the same 
object (Figure 5D). This latter category serves as one possible definition of a 'shared cause' for the two edges. 

A second definition of shared cause invokes not just the object membership but also the relative depth of the 
objects. In particular, we use ordered membership functions □ (Section 4.5), and wc classify these □ according to 
whether a pair of pixels from each edge both falls into the same object and that object is above the object present 
at the remaining pixels (Figure 5E). The relevant □ are listed in Table IB. 

With either definition, the likelihood ratio of shared cause to different cause is _L = X^aes -^('^)/X]aeD ^('^)' 
where S and D are the sets of membership functions categorized as shared or different causes respectively. This 
likelihood ratio varies as a function of the positions and relative orientation of the two edge pairs (Figs. 5C-D). 



4.5 Generalizations 

We can calculate the relative depth of objects by using an ordered object membership function □ rather than an 
unordered membership function n. (The Hebrew letter final mem D is used only at the end of a word, representing 
that object order matters.) Whereas Q was a set of subsets, □ is an ordered set of subsets with □„ representing the 
pixels contained by the nth-highest object sampled by any of the N selected pixels. The recursion in this case is 
even simpler than Equation 3: 

J- — ^a(u,a) 

There is no summation here because there is only one term for which the first object is highest in the stack of 
objects. One may use a partial ordering if not all relative depths are of interest, and then there will be a sum over 
arrangements consistent with the partial ordering. 

Note that there are more hidden variables of interest besides the object membership and relative depth, and the 
joint probabilities of these can be calculated by a similar recursive formula, without marginalizing away the hidden 
variables. The joint distribution of shape and membership, for instance, can be calculated as 



P{^,P) = 1 n y2PiPn)Qcr(n,n){pn)P{J^\n,P\n) 



where p is now a vector of N shape parameters, with pn indicating the shape parameters for the topmost object 
present at pixel location x„. 
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A: Classification of unordered a 



S: Shared cause 


D: Different causes 


1256 1 34 1 78 
1278 34 56 
12 1 56 1 3478 
12 78 3456 
1256 3478 
1278 3456 


12 1 34 1 56 1 78 



B: Classification of ordered □ 



S: Shared cause 


D: Different causes 


1256 > 34 1 78 
3478 > 12 56 
1278 > 34 56 
3456 > 12 78 
1256 1 3478 
1278 3456 


34 1 78 > 1256 34 > 1256 > 78 78 > 1256 > 34 
12 56 > 3478 12 > 3478 > 56 56 > 3478 > 12 

34 56 > 1278 34 > 1278 > 56 56 > 1278 > 34 
12 78 > 3456 12 > 3456 > 78 78 > 3456 > 12 
12 1 34 1 56 1 78 



Table 1: Object membership functions used for joint edge statistics. For compactness wc represent object mem- 
bership functions by the pixel indices divided symbolically into ordered or unordcrcid groups. For example, 
{{xi,X2}, {x3,X4}} is written as 12|34 if unordered, and as 12 > 34 if ordered such that the object containing 
points xi and X2 lies above the object containing X3 and X4. These object membership functions are classified ac- 
cording to whether they reflect a shared cause or different causes for the two edges, using unordered (A) or ordered 
(B) representations. 

4.6 Empirical sampling of dead leaves and natural images 

For probabilities involving many image points, we generate many dead leaves images and empirically sample from 
them to obtain histograms. Images are produced by layering objects from front to back until all image pixels are 
members of some object, a process that yields stationary image statistics [46]. 

Natural images were drawn from van Hateren's image database [47]. Feature distributions were obtained by 
log-transforming images [12], filtering them by the relevant Haar wavelets, and computing univariate or bivariate 
histograms. 
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Figure SI: Joint probabilities of pixel intensities, based on an ensemble of elliptical objects and gaussian-distributed 
object intensities with an additive gaussian white noise texture (Methods). Contour plots are shown for two pixels 
(A) and three pixels arranged in an equilateral triangle (B). These joint distributions are weighted averages of 
independent and correlated distributions. The weighting factors are the various object membership probabilities 
which are plotted below the joint intensity distributions as a function of the distance between pixels. 
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Figure S2: All object membership probabilities -P(!3) as a function of the spacing between the pixels, for four 
pixels arranged in a square (right). Two different shape ensembles are shown, circles and ellipses, with sizes given by 
P{r) oc for r ranging between the limits r_ and r+ (dashed lines). Curves are labeled by their object membership 
functions. Symmetrically permuted membership functions have identical curves. For circles, the configuration {14|23} 
is impossible, but otherwise the curves for circles and ellipses are remarkably similar across all pixel spacings, because 
both shape ensembles have similar local properties (extended edges) and global structure (convex shapes with the 
same size distribution). Since objects have sharp edges that closely spaced pixels rarely straddle, nearby pixels almost 
always fall into the same object, with P({1234}) w 1. When pixel spacing exceeds the largest object dimension, 
no two pixels can fall into the same object, so the only membership function allowed is a = {1|2|3|4}. With pixel 
spacings between these extremes, many more object membership probabilities take on nonzero values. 
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Figure S3: Mixture distributions and mixture components of sparsely sampled Haar wavelet features, calculated 
exactly for dead leaves images. (A,B) Contours of the log-probabilities logP(/i,/2) for colocalized, orthogonal 
wavelets /i and /2, sampled with four or eight points per feature (insets). (C.D) Elliptical contours of jointly 
gaussian mixture components P(/i,/2|!2), shaded according to their weight The mixture distributions already 

have rounded diamond contours formed from weakly correlated components, as well as some strongly correlated and 
anti-correlated components which appear at all angles with dense sampling (Figure 3B) . (E-H) Joint log-probabilities 
and mixture components for nearby, parallel Haar wavelets, plotted as in A-D. The anticorrelation and 'pillow' shape 
of these distributions are already visible with sparse sampling of the features. 




Figure S4: Detailed geometry for Figure 5. (A) An edge exists when an object splits four pixels into two pairs. 
Pixels arranged in a rectangle with an aspect ratio of Ax/ Ay = 2.75 permit a range of edges with a 40° orientation 
bandwidth as used in [20]. (B) Pairs of edges thus defined are related by three parameters: distance d, orientation 
difference 9, and relative direction 0. 
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Figure S5: Diagrams for illustrating inclusion probabilities Qa-{a,n)- A 'leaf function showing the shape of an 
object. Li(x,p) = 1 at points x that are inside an object of shape parameter p, and Lq(:x., p) = 1 at points outside 
it. Here the shape parameter p specifies a smooth irregular object. (B) Example indicator functions Q<T(B,n) (c, p) 
identify locations c where an object could be placed to enclose all pixels £ a„ and exclude the rest. Rotated copies 
of the object shape surround each point x^, designating the locations c where that object will enclose x^. An arrow 
points to the shaded region where an object could be placed to enclose both xi and X3 but not X2, whose area is 
<9ioi(p)- The other shaded region indicates locations where an object would enclose only X2, whose area is Qoioip)- 
Note that this diagram represents possible locations c of a single object, not three objects! 




Figure S6: Simplified representations of object configurations that dominate feature distributions at high amplitudes. 
Colors indicate objects of unspecified intensity, shading indicates weighting by Haar wavelets. (A) A Haar wavelet 
Hi takes a difference of intensities Ii^i and l2,i each totalled over a finite region. The pixels Q^'* contained both 
in elliptical object a„ and in region 2 of wavelet i are outlined. (B) Colocalized, orthogonal Haar wavelets with 
circular support. (C, D) Parallel, nearby Haar wavelets, with objects that induce negative and positive correlations, 
respectively. To simplify the calculations, objects differ only in their horizontal extent, and extend completely to 
either the left or right edge of each wavelet. The relevant variable is then the width of the overlap between the object 
and the wavelet filter, denoted di and dr. Compare these simplified configurations to those shown in Figures 3C and 
3G,H. 
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SI Inclusion probabilities for an ensemble of ellipses 



In the main text we reported a universal recursion equation expressing object membership probabiUties -P(!3) in terms 

of some geometric factors Qa-{a.n) which depend on the shape ensemble. There we showed how these probabiUties 
could be expressed geometrically, by first averaging Qa{n,n){^j p) over position c via contour integrals, and then 
averaging over the shape ensemble p. Here we explain in detail how inclusion probabilities Q„(a,n) can be calculated 
exactly for an ensemble of circular objects. We then use a simple transformation to generalize the result for circles 
to an ensemble of ellipses. When the dust settles, we will have averaged Qcr{c,p) over positions c and shapes p and 
for all binary vectors <r. 

For an ensemble of circles, the shape parameter p is jiist a radius r, which we draw ivom a scale-invariant size 
distribution P{r) oc r~^. Circular contours are easy to express analytically. However, as described in the main text, 
the integrals of Q„{c,r) over both the contours and size ensemble are more difficult because they must be done 
piecewise. We do this in two steps. First, we evaluate the general form of the indefinite integrals at the cndpoints 
of the piecewise intervals. Second, we describe an algorithm that synthesizes these isolated contributions into the 
complete piecewise integral, yielding the desired Qcr- 



Sl.l Parameterizing circular contours 

Equation 4 related the positional average Qa{a.n){p) to the total area of the region where QCT(D.n)(c,p) = 1, and 
thence to a contour integral. In this section we evaluate this contour integral for circles with fixed radius, so that 
p = r. It is helpful to change from the generic notation used in Section 4.1 to a notation which is specific to circular 
objects. As shown in Figure S7A, the boundaries of regions with constant (5tT(a,n)(c, r) are all circular arcs centered 
on some point Xj, 

Si{t) = ret + s^i 

with a unit vector defined as et = (cost, sin f). Each arc terminates at angles t of the form 

tij± = Oij ± (pij = tan~^ (xj -Xj)± cos~^ {uij/'2r) 

where 9ij is the angle of the line connecting the circle centers, and ±(pij are the angles that the intersection points 
make with that line (Figure S7A). 9ij is independent of r, whereas depends on the ratio of r to the distance 
Uij = \xi — Xjl between the circles as = cos~^ {uij/2r). 



SI. 2 Contour integration 

Since the unit normal vectors are simply h{t) = et and the arc length is ds = |s(t)|(if = rdt, we can now easily 
perform the contour integral (Equation 5) over each arc analytically. 

Am = 7y I ^m{t) ■ n(t) ds = - I (r^ + rxi • e*) dt = aik±{r) - aij±(r) 



where we have defined 



rOij r^ _^Uij Uij ^ r , ^ij 

= „ ± — cos — ^ H — T^Xj • ee ._iL ± -Xi • ee-.l/ 1 - — ^ (SI) 
2 2 2r 4 2 2 V ^ ^ 

For r smaller than the distances between pixels, the circular arcs do not intersect and are thus complete circles with 
total area of nr^, as expected. 

Finally we can obtain the total area (3o-(d,ti)(^) by adding up the relevant aij±{r) appropriately. We defer 
discussion of this step to Section SI. 5. 



SI. 3 Indefinite integral over radius 

Next we have to average these quantities over the distribution of r. To achieve scale invariance in generated images, 
the distribution of object radii P{r) should be proportional to [12]. Some deviation from this scaling behavior 
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Figure S7: (A) Diagram depicting the quantities needed to calculate inclusion probabilities Q a- {a, n) {''')■ The different 
regions of constant Qa-{a,n){c,r) for fixed r are bounded by circular arcs centered on the pixels x. Highlighted is one 
particular arc Si{t) centered on point x^. This arc is bounded by tij- and ti^^, two angles at which other circles 
intersect. Centers and x^ are separated by the distance Uij and angle 9ij. The location at which the corresponding 
circles intersect deviates from the line connecting the centers by angle —4>ij, so that Uj- = 9ij — ipij. (B) Illustration 
of how the contours around regions with constant Qrr(a,n){Cif) change shape as r increases (from left to right). Two 
regions in c-space first touch when r equals half the distance between two pixels x^ and Xj , a critical radius r*j we 
call a 'kissing point' (center panel). As r increases further a new contour of reversed orientation is created, bounding 
a region within which an object of radius r can enclose both pixels. (C) Similarly, a 'triple intersection' always 
exists for a particular r*jj,, the circumradius, at which any three non-coUinear pixels x^, Xj and x^ are equidistant 
from a fourth point called the circumcenter (center panel). As r crosses this critical radius, the existing contour 
connecting the three intersection points changes orientation, and the enclosed region is associated with a different 
QrT(-a,n)- (D,E) Illustrations of two strategies for integrating Qcr(r) over r: choose only one region at a time, and 
track its contour as r varies (D); or track all contour endpoints over r and add their contributions to all appropriate 
regions (E) . We use the latter strategy. Arrows in panel F depict the four regions cr receiving identical contributions 
(up to a sign) from the contours along Si{t) that terminate at the intersection point xi3_. 
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is required to prevent images from degenerating with high probabihty into white noise or uniform coloring [12, 28]. 
Here we choose to set upper and lower size cutoffs r e [r_,r+] to satisfy this constraint, so that 



P(r) 



3 J, g [r_,r+] 
otherwise 



with z = \{rZ^ - r+^). 

The r-dependence of aij±(r) in Equation SI takes the forms r\j 1 — {u/2r)'^ , r^cos"^ {u/2r), and r^. Each of 
these terms must be averaged over P{r). The first average can be solved analytically as 



/ 



P{r)rJl- {u/2rfdr =- [ r^^Jl- {u/2rf = -- 



1 • ^ , 7^ 
-sm2<^+- 



The average of the second term can be found in tables of integrals, and involves a special function, the dilogarithm 
Lhiz). 

j P{ry cos-1 ^dr=^^jr-^ cos"! |; dr = ^</)2 - log (l + e^'^) + Lis (-e''^) 

For u < 2r (required for the two relevant circles to intersect), the imaginary component is constant and therefore 
cancels in any real definite integral. We can therefore take just the real component without influencing the result. 



/ 



P{ry cos-^ ^dr = -<^log ^ - ^^[iU2 (-e^*^)] 



The imaginary part of the dilogarithm evaluated on the complex unit circle is related to another special function 
known as Clausen's integral, for which optimized numerical routines have been written [48]. 

^[iU2 (-e^'"^)] = CI2 (-2<^ - tt) 



The remaining terms in aij±{r) are elementary to integrate: / P{r)r^dr = | logr and / P{r)dr — 

Combining all these pieces with their correct coefficients, we obtain the indefinite integral for the size average of 
aij±{r). 



/I u 
dr P{r)aij±{r) =- —6' logr + ^^Xi-eg. 



SI. 4 Identifying piecewise smooth intervals over radius 

The definite integral over r must be performed piecewise because its integration contours may change at certain 
critical radii r*. Generically, there are two types of critical radii, depicted in Figure S7B,C: 'kissing points' where 
r*j is half the distance Uij between a pair of points and xj, so that their corresponding circles just touch; and 
'triple intersections' where r*^^, equals the cireumradius of three points x^, Xj and x^., so that the three corresponding 
circles all meet. For three points separated by distances Uij, ujk, and Uki, and semiperimeter s = -^{uij +Ujk +Uki), 
the cireumradius is 

'^ijk = UijUjkUki/4:^ S{S - Uij){s - Ujk){s - Uki) 

If the pixel locations have extra symmetries, e.g. lie on a lattice, then several critical radii r* may coincide. In 

this case each r* can be treated sequentially without changing the result, as if perturbing each r* infinitesimally: 
t'ij±i'>'") ~ bij±{r') contributes zero in the limit r" — r' — )■ when there are no intervening critical radii. 
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SI. 5 Mapping piecewise integrals onto appropriate Qcr 



Now we must calculate Qa- by adding up the definite integral bij±{r) evaluated at the appropriate critical radii r* 
and the relevant triples {i,j,±). Consider two strategies for this. First, one could choose one particular cr, and 
track how the cusps of Q„-{c,rys boundary appear, change, and disappear as a function of r, and then add up the 
appropriate contributions from Equation S2 (Figure S7D). One would then repeat this procedure for every possible 
cr. Second, one could choose a particular intersection point Xjj-i- between two objects, track how it is associated with 
different regions as a function of r, and add its contribution to the various appropriate By iterating through all 
intersection points, eventually all contributions to all Qrr arc computed (Figure S7E). This latter strategy is easier 
because the behavior of the intersection points is simpler to track than the various (possibly unconnected) regions 
where Q^{c,r) = 1. This is the approach we describe below. 

To compute the definite integral corresponding to Eqiiation S2 above, we must therefore associate each integrand 
aij±{i') with boolean vectors cr designating the correct targets Q„ for each interval of r. The region geometry, and 
thus these desired associations, change only at critical radii; between critical radii the associations are constant. By 
construction, aij±{r) (Equation SI) is the result of a contour integral terminating at an intersection between circles 
centered on Xj and x^- (Figure S7A). We label this intersection point by Xij± = x, + re0..±^... Contour integrals 
terminating at this point contribute to every one of the four regions that touch Xij±, i.e. the a involving all four 
allowed combinations of its elements (Tj e {0, 1} and <jj € {0, 1} (Figure S7F). The point Xjj± is not on the boundary 
of any circles centered on other pixels x^, since otherwise there would be a critical radius within the selected r 
interval. Xjj± is thus either strictly inside or strictly outside a circle of radius r for all f ^ i, j. We can now specify 
all elements of cr as 



where ii(xjj± — x^,r) is the leaf function from Methods Section 4.1. This relation identifies the appropriate targets 
Qrr for the hij± of Equation S2. 

To identify the signs Cij± with which bij± contribute to the target Qcr, it helps to go back and compute the 
signs that aij±{r) contribute to the target area Qtrir). These signs depend on the geometry of the region contours. 
Consider how the region boundaries change their geometry as r increases from r_ to r_|_. A contour around an object 
boundary is counterclockwise initially, i.e. before the contour intersects any other object boundaries. As r increases 
past a kissing point r*j , a pair of intersection points Xjj± is created along with a new region with clockwise orientation 
(Figure S7B). Note that the contours around the object centered on x^ initially converge at an intersection x,j_ and 
diverge at Xjj_|_. In other words, intersections Xjj-_ are initially endpoints of the contours along Si{t) that contribute 
+aij± to the contour integral (Equation 5), and Xij+ are initially starting points that contribute —aij±. However, 
as r increases past each triple-intersection r*^^, for k ^ i,j, another circle centered on x^ encloses the intersection 
point. The orientations of the contours at Xij± then reverse (Figure S7C), and the sign that each aij± contributes 
also reverses. Thus the overall convergence for paths at an intersection point is: converging for — , diverging for +, 
and reversed by the number of circles enclosing the point. Mathematically, we can write the desired sign as 



Note that Cij±{r) does not vary between critical radii r*, so we may use its value anywhere within the integration 
interval. Finally, when we integrate aij±{r) over r' < r < r", the value of the indefinite integral bij± at r' is 
subtracted from the value at r". Thus, for each interval between critical radii we add 



to the appropriate Qa-- 

There is one remaining subtlety in adding up the contributions to Qa-- In the first term of bij± there is an 
ambiguity of 2n in what angle is subtended by a given arc, which cannot be resolved by local properties of the arc 
endpoints alone. We remedy this by computing AQa{r' ,r") modulo ^ log r"/r', which is the maximum possible 
contribution an area can make between r' and r" . This guarantees that we update Qo- with the unique definite 
integral over r' < r < r" that lies between and this maximum. 

SI. 6 Sumniciry of the algorithm for calculating Q^- 

This completes the mathematics necessary to calculate the Qcr- To summarize, we present the method in algorithmic 
form. 




Cij±{r) = ^(-l)i:.^^.,i'i(-^-x«±,r) 



AO.(r' 



•',r") = Q,±(^)-(6,,±(r")-5,,± 



(r')) 
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1. Initialize all Qg^ to zero. 

2. Add JJ^* drP{r)nr'^ = ^ log ^ to Qg^ for each circle, where r* = m.mj^ir*j is the first kissing point for that 
circle and Si is a vector of zeros with a 1 at index i. This is the area accumulated in Q^. before any other 

circles were touched. 

3. Sort all critical radii r*j and r*^^, within the integration bounds r_ and r+. 

4. For each interval r' < r < r" bounded by sequential critical radii: 

(a) For each existing intersection point Xij±: 

i. Calculate the region indicators cr to which the point Xjj± contributes 

ii. Add AQa.{r',r") modulo f log ^ to Q^r 

5. SetQo = l-E^^o<3-- 

Once the Qcr(ii,n') a-J't- calcidatcd for all object membership functions a, then Qa-{n.^,^.k) must be calculated for the 
reduced used in the recursion. For efficiency, this can be accomplished by marginalizing Qo- over the appropriate 
indices <Ji, rather than recalculating it with a smaller set of pixels. 

Note that with a different size ensemble P{r), the expression for bij± would change, but the procedure for 
combining them to obtain the Qcr would be the same. 

SI. 7 Converting from circles to ellipses 

It is straightforward to transform our calculation of Qa- for circles into a result for ellipses of equal area but eccentricity 
e and orientation tp. All distances are effectively scaled by y/e in the direction of and l/\/e in the orthogonal 
direction. This is equivalent to transforming the pixel locations 

1 / e cos^ + sin^ -0 (e — 1) cos V'sin-0 
* \ (e — 1) cos sin^/) cos^ V + e sin^ 

and recomputing the Qcr with these x^(e,7/')- 

Unfortunately we cannot analytically integrate bij±{r*) as a function of eccentricity e or angle ip, because the 
dependence on the points x.^{€,■tp) already involves special functions. Instead, to obtain the average over possible 
ellipses we use a discrete ensemble of eccentricities and angles and sum over them as 

(Q^(c,r,e, V))c,.,e,^ = 5^Q,(e,V)P(e)P(V) 

More generally, when the integral cannot be expressed analytically using easily computable functions, one may 
specify the ensemble P{p) by a discrete number of allowed shapes, and compute the ensemble average as a sum 
rather than as an integral. 

The result of these calculations are concrete numbers for the inclusion probabilities Qcr, which can then be 
substituted into Equations 1, 2, and 3 to calculate the object membership probabilities and joint distributions of 
pixel intensities and image features. 

References for Supporting Information 
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Hebrew letter mem: an object membersmp lunction 


|q| 


number oi distinct objects m a 


On 


set of pixels contained in nth object 




object membership function with rath object removed 


Pa 


probability that pixels are divided according to Q 


X 


locations of all iV selected pixels 


Xj 


location of ith pixel 


1 


vector of all N pixel values 


T 


pixel value at point Xj 


T 


vector of all pixel values in nth object 


crip, n) 


boolean vector indicating pixels in nth object 


W<7{0,0) 


probability that an isolated object includes no selected pixels 


/-") 

Vcr(a,n) 


probability that an isolated object includes only pixels 0^ 


V(T(a,rt) 


as above, but given the object shape p 




as above, but also given the object position c 




location of two intersections between objects centered on Xj and Xj 




Qisxance oeLween Xj ana x^ 


ft ■ 


dllgiC Ui Lilt; VcCLUl Ji.j — Ji-i 


A- ■ 
(pij 


absolute value of angle between intersection points and line connecting and Xj 








unit vector (cos 0, sin 0) 




contour around object centered on Xj 


aij±{r) 


indefinite integral over t along contour Sj(t) evaluated at tij± with fixed r 


bij±{r) 


indefinite integral of aij±{r) over r 


Cij±{r) 


sign indicating whether contour Si{t) starts or ends at t = tij± 




critical radius at which regions of constant Qa-{a,n) (c, p) change structure 


ij 


radius of kissing point for circles on x, and Xj 


r* 
ijk 


radius of triple intersection for circles on Xj, x^ and 


e 


ellipse eccentricity 




orientation of major axis of ellipse 




transformed pixel location 



Table SI: Glossary of symbols used 
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